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Abstract. The present study is devoted to the problem of tsunami wave generation. The 
main goal of this work is two-fold. First of all, we propose a simple and computationally 
inexpensive model for the description of the sea bed displacement during an underwater 
earthquake, based on the finite fault solution for the slip distribution under some assump- 
tions on the dynamics of the rupturing process. Once the bottom motion is reconstructed, 
we study waves induced on the free surface of the ocean. For this purpose we consider three 
different models approximating the Euler equations of the water wave theory. Namely, we 
use the linearized Euler equations (we are in fact solving the Cauchy-Poisson problem), 
a Boussinesq system and a novel weakly nonlinear model. An intercomparison of these 
approaches is performed. The developments of the present study are illustrated on the 
17 July 2006 Java event, where an underwater earthquake of magnitude 7.7 generated a 
tsunami that inundated the southern coast of Java. 



1. Introduction 

Tsunami waves have attracted a lot of attention by researchers. The interest of the 
scientific community has especially increased after the two megatsunamis in December 2004 
[SB06], where nearly 230,000 people in fourteen countries lost their lives, and in March 
2011, where 20,000 people lost their lives in Japan. The 2004 event also led Indian Ocean 
countries to develop Tsunami Warning Systems (TWS) [Syn05, Bas06], unfortunately more 
on an individual basis than on a collective basis. The most elaborated warning system to 
date is the Pacific Ocean TWS, which has been developed over several decades by efforts 
of NOAA's specialists [TGB+05, GBM+05]. 

An operational tsunami wave modeling tool is an essential part of any warning system 
[TGB + 05, TDS07]. Mathematical and numerical models in use should be constantly im- 
proved to produce more accurate results in less CPU time [Ima96, TG97, DPD11]. In order 
to study the propagation of a tsunami wave, an initial condition must usually be provided 
to any numerical model designed for this purpose. The present study is an attempt to 
improve the construction of the initial tsunami waveform. The set of existing practices de- 
scribed in the literature constitutes the field of the so-called tsunami generation modeling 
[Ham73, TT01, DD07d, Dut07, DD09, DD10]. 

The modeling of tsunami generation was initiated in the early sixties by the prominent 
work of Kajiura [Kaj63], who proposed the translation of the static sea bed displacement 
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towards the free surface as an initial condition. Classically, the celebrated Okada [Oka85, 
Oka92] and sometimes Mansinha & Smylie 1 [MS67, MS71] solutions are used to compute 
the co-seismic sea bed displacements. This approach is still widely used by the tsunami 
wave modeling community. However, significant progress has been recently made in this 
direction [OTM01, DD07d, Dut07, DD09, RLF+08, SF09, DPD11]. 

In the present study we exploit some recent advances in seismology to reconstruct better 
co-seismic displacements of a tsunamigenic earthquake. More precisely, we suggest using 
the so-called finite fault solution developed by Ji and his collaborators [BLMOO, JWH02], 
based on static and seismic data inversion. This solution provides multiple fault segments 
of variable local slip, rake angle and several other parameters. By applying Okada's solution 
to each subfault, we reconstruct the sea bed displacement with higher resolution. To our 
knowledge, this technique has already been employed to model the Kuril islands tsunamis 
of 15 November 2006 and 13 January 2007, cf. [RLF^OS]. Since Okada's solution consists of 
relatively simple closed-form analytical expressions, all computations can be done efficiently 
enough so that they can be used in a real-time TWS (cf. [WL08]). The obvious sine qua 
non condition is that the corresponding finite fault inversion should also be performed in 
a reasonable time. 

In the present study we go further in reconstructing the dynamic sea bed displacement 
according to the rupture propagation speed and the rise time also provided by the finite 
fault solution. Constructed in this special way, sea bed displacements are then coupled with 
several water wave models. Among them, there is a novel weakly nonlinear solver based on a 
formulation involving the Dirichlet-to- Neumann operator which is computed approximately 
using Fourier transforms. The other two models considered here are the linearized free 
surface Euler equations and a Boussinesq type system. Developments presented in this 
paper are illustrated on the example of July 17, 2006 Java event [AKLV06]. However, 
we would like to stress that the methodology presented in this study is quite general and 
can be applied to many other tsunamigenic earthquakes for which a finite fault solution is 
available. 

The paper is organized as follows. In Section 2 we describe the static and dynamic sea 
bed displacements, while in Section 3 we present a simple approximate water wave solver 
with a moving bottom. In Section 4 we study numerically the generation process of a 
real- world event. An intercomparison of the three models mentioned above is performed. 
Some important conclusions are drawn in Section 5. 

2. CO-SEISMIC DISPLACEMENT CONSTRUCTION 

The modeling of tsunami generation is directly related to the problem of the bottom 
motion during an underwater earthquake. Traditionally, Okada's solution [Oka85, Oka92] is 
used in regimes characterized by an active fault of small or intermediate size, i.e. consisting 
of one or a few segments (e.g. the great Sumatra 2004 earthquake, [SB06, IAK + 07]). In 
this case the resulting vertical displacement field is translated to the free surface. This 
approach is conventionally referred to as passive tsunami generation [DDK06], contrary to 



In fact, the Mansinha & Smylie solution is a particular case of the more general Okada solution. 
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Fault width, km 


40.0 


Focal depth, km 


20.0 


Slip, m 


2.5 


Dip angle 


10° 


Slip angle 


95° 


Strike angle (clockwise from N) 


289° 



Table 1. Seismic fault parameters for the Java 2006 event. The correspond- 
ing seismic moment can be taken as Mq = 2.52 x 10 27 N- m (M w = 7.56). 

the active generation which explicitly involves the bottom motion dynamics [DD07d]. Since 
our methods will be illustrated on the example of the July 17, 2006 Java event, we show in 
Figure 1 a typical single-fault based initial condition used for the corresponding tsunami 
wave modeling [Yal08]. The seismic parameters used to produce this vertical displacement 
are given in Table 1. 

Remark 1. The celebrated Okada solution [Oka85, Oka92] is based on two main ingredients 
- the dislocation theory of Volterra [Vol07] and Mindlin's fundamental solution for an 
elastic half-space [Min36]. Particular cases of this solution were known before Okada's 
work, for example the well-known Mansinha & Smylie's solution [MS67, MS71]. Usually, all 
these particular cases differ by the choice of the dislocation and Burger's vector orientation 
[Pre65] . We recall the basic assumptions behind this solution: 

• The fault is immersed into a linear homogeneous and isotropic half-space 

• The fault is a Volterra type dislocation 

• The dislocation has a rectangular shape 

For more information on Okada 's solution we refer to [DD07d, DD07a, Dut07]. 

The finite fault solution is based on the multi-fault representation of the rupture [BLM00, 
JWH02]. The rupture complexity is reconstructed using a joint inversion of the static and 
seismic data. The fault's surface is parametrized by multiple segments with variable local 
slip, rake angle, rise time and rupture velocity. The inversion is performed in an appropriate 
wavelet transform space. The objective function is a weighted sum of L±, L 2 norms and 
some correlative functions. With this approach seismologists are able to recover rupture 
slip details [BLM00, JWH02]. This available seismic information is exploited in this study 
to compute the sea bed displacements produced by an underwater earthquake with higher 
geophysical resolution. 

The proposed approach will be directly illustrated on the Java 2006 event. The July 17, 
2006 Java earthquake involved thrust faulting in the Java trench and generated a tsunami 
wave that inundated the southern coast of Java [AKLV06, FKM + 07]. The estimates of 
the size of the earthquake (cf. [AKLV06]) indicate a seismic moment of 6.7 x 10 20 N- 
m, which corresponds to the magnitude M w = 7.8. Later this estimate was refined to 
M w = 7.7 [Ji06]. Like other events in this region, this 2006 event had an unusually low 
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Figure 1. Static vertical displacement in meters of the seabed computed 
with the single fault parameters provided in Table 1. The maximum lift is 
0.7215 m while the maximum subsidence is 0.4030 m. The x— axis is the 
longitude while the y— axis is the latitude. The y— axis points to the North. 

rupture speed of 1.0 - 1.5 km/s, and occurred near the up-dip edge of the subduction 
zone thrust fault. According to Ammon et al, most aftershocks involved normal faulting 
[AKLV06]. The rupture propagated approximately 200 km along the trench with an overall 
duration of approximately 185 s. The fault's surface projection along with ocean ETOPOl 
bathymetric map are shown in Figure 2. We note that the Indian Ocean bathymetry 
considered in this study varies between 7186 and 20 meters in the shallowest region. 

Remark 2. The estimate of the finite fault inversion for this earthquake was also performed 
by the Caltech team [Ozg06]. The magnitude estimated in that study was M w = 7.9. In this 
study we do not present numerical simulations using their data but it is straightforward to 
apply our algorithms to this case as well. 

2.1. Static displacement. In order to illustrate the advantages of the proposed approach 
we will also compute the static co-seismic displacements using the finite fault solution [Ji06]. 
The fault is considered to be the rectangle with vertices located at (109.20508° (Lon), 
-10.37387° (Lat), 6.24795 km (Depth)), (106.50434°, -9.45925°, 6.24795 km), (106.72382°, 
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(a) Top view. (b) Bathymetry side view. 

Figure 2. Surface projection of the fault's plane and the ETOPOl bathy- 
metric map of the region under investigation. The symbol * indicates the 
epicenter's location at (107.345°, —9.295°). The local Cartesian coordinate 
system is centered at the point (108°, —10°). The region is located between 
(106°, -8°) and (110°, -12°). 



P-wave celerity c p , m/s 



S- wave celerity c s , m/s 



Crust density p, kg/m 3 



Dip angle, 5 



Strike angle (clockwise from N) 
Table 2. Geophysical parameters used to mode 
subduction zone in the region of Java. 



6000 



3400 



2700 



10.35° 



289° 

elastic properties of the 



-8.82807°, 19.79951 km), (109.42455°, -9.74269°, 19.79951 km) (see Figure 2a). The 
fault's plane is conventionally divided into N x = 21 subfaults along strike and N y — 7 
subfaults down the dip angle, leading to a total number of N x x N y = 147 equal segments. 
Parameters such as subfault location (x c ,y c ), depth di, slip u and rake angle </> for each 
segment are given in Table 3 and can also be downloaded from [Ji06]. The elastic constants 
common to all subfaults and parameters such as dip and slip angles are given in Table 2. 
(We note that the slip angle is measured conventionally in the counter-clockwise direction 
from the North. The relations between the elastic wave celerities c p , c s and the Lame 
coefficients A, /i used in Okada's solution are given in Appendix C.) 

We compute Okada's solution at the sea bottom by substituting z = in the geo- 
physical coordinate system and taking the vertical component of the displacement field 
Oi(x; 5, A, /i, . . .), where 6 is the dip angle, A, \i are the Lame coefficients (see Appendix C) 
and the dots denote the dependence of the function 0(x) on eight other parameters, cf. 
[DD07d]. The resulting co-seismic vertical bottom displacement can be computed as 
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Figure 3. The vertical displacement of the finite fault solution, cf. [Ji06]. 
The corresponding seismic moment is M = 3.53 x 10 27 N- m (M w = 7.65). 
The maximum lift is 0.4629 while the maximum subsidence is 0.1997. 

a simple superposition of subfault contributions: 

N x xN y 

i=i 

The graph of ((x) is presented in Figure 3. The specific static displacement can be com- 
pared with the single fault classical approach depicted on Figure 1. It is worth mentioning 
that more than one local extrema can be found in this solution due to a higher slip reso- 
lution. 

Hereafter we will adopt the short-hand notation Oi{x) for the vertical displacement 
component of Okada's solution for the i th segment having in mind its dependence on 
various parameters from Tables 2 and 3. 

2.2. Dynamic co-seismic displacements. Here, we go even further in the reconstruc- 
tion of the bottom motion. By making some assumptions on the time dependence of the 
displacement fields, we can have an insight into the dynamics of the sea bed motion. 
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The finite fault solution provides two additional parameters concerning the rupture dy- 
namics for the July 17, 2006 event — the rupture velocity v r = 1.1 km/s and the rise time 
t r = 8 s. The epicenter is located at the point x e = (107.345°, —9.295°) [Ji06]. Given the 
origin x e , the rupture velocity v r and the i th subfault location Xi (the full list is provided 
in Table 3), we define the subfault activation times t{ needed for the rupture to reach the 
corresponding segment % by the formula: 

U = — — — , i = 1,...,N X x N y . 

For the sake of simplicity and due to the lack of information we assume implicitly that the 
rupture speed v r is constant along the fault; however this can be refined in future studies. 

We will also follow the pioneering idea of J. Hammack [Ham72, Ham73] developed later 
in [TT01, THT02, DD07d, DDK06, KDD07] where the maximum bottom deformation is 
achieved during some finite time (known as the rise time) according to a specific (in an ad 
hoc manner) dynamic scenario. Various scenarios on the time dependence (instantaneous, 
linear, trigonometric, exponential, etc) can be found in [Ham73, DDK06, DD07d]. In this 
study we will adopt the trigonometric scenario which can be described by the formula: 

T(t) = H(t - t r ) + ]-U{t)U{t r - t) (l - cos(vrf/t r )) , 

where %{t) is the Heaviside step function. For illustrative purposes this dynamic scenario 
is represented on Figure 4. Physically the function T(t) represents the time history of the 
vertical bottom displacement in terms of its final amplitude. We assume that during the 
rise time temporal interval [0,t r ] the vertical displacement goes from zero to its final stage 
according to the trigonometric scenario. 

Finally, we put together all the ingredients in order to construct the dynamic sea bed 
motion: 

N x xN y 

C(x,t)= nit-t^Tit-t^Oiix). (2.1) 

1=1 

In the following sections we will present several approaches to couple this dynamic de- 
formation with the hydrodynamic problem to predict waves induced on the ocean's free 
surface. 



3. Fluid layer solution 

Once the sea bed deformation is determined, a water wave problem must be solved in or- 
der to compute the free surface motion induced by the ocean bottom shaking. Traditionally 
this difficulty is circumvented by the simple translation of the static bottom deformation 
onto the free surface [Kaj63], known as the passive generation approach [DDK06, KDD07]. 
In this section we present three approximate models to the water wave problem with mov- 
ing bottom that we will use in combination with the finite-fault solution to study the 
tsunami generation problem. 
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Trigonometric scenario 




Figure 4. Trigonometric scenario with rise time t r — 1 s. 



3.1. Linearized Euler equations — CP model. Consider an ideal incompressible fluid 
of constant density p. The horizontal projection of the fluid domain Q is a subset of 
M. 2 . The horizontal independent variables are denoted by x = (x,y) and the vertical 
one by z. The origin of the cartesian coordinate system is chosen such that the surface 
z = corresponds to the still water level. The fluid is bounded below by the bottom 
z = —h(x,t) and above by the free surface z = r](x,t). Usually we assume that the total 
depth H(x,t) := h(x,t) + rj(x,t) remains positive H(x,t) > h > at all times t G [0,T]. 
The sketch of the physical domain is shown in Figure 5. 

Remark 3. Classically in water wave modeling, we make the assumption that the free 
surface is a graph z = T](x,t) of a single-valued function. It means in practice that we 
exclude some interesting phenomena, (e.g. wave breaking phenomena) which are out of the 
scope of this modeling paradigm. 

The linearized water wave problem consists of the following set of equations [Ham 72, 
Ham73, DD07d]: 



dtV 
d t (f> + gr) 
d t h + d z (f) 



0, (x,z) e n x [-h,0], 

0, z = 0, 

0, z = 0, 

0, z — —h(x, t). 



(3.1) 
(3.2) 
(3.3) 
(3.4) 



This set of equations together with an initial condition is also often referred to in the liter- 
ature as the Cauchy-Poisson (CP) problem after the pioneering work of Cauchy [Cau27]. 
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Figure 5. Sketch of the physical domain. 



In view of the specific requirements of the analytical techniques used in the applications, 
we will assume first that the domain Q = IR 2 , i.e. it is unbounded in the horizontal extent, 
and the bottom has a special form: 

h(x,t) =h - ((x,t), 

where h is some uniform depth and ((x, t) is the sea bed displacement due to an underwater 
earthquake. In Section 2.2 one possible construction of the bottom displacement was 
proposed. Using integral transform methods (cf. [Ham 73, TT01, DD07d, KDD07]), one 
can derive the following expression for the free surface elevation n(x,t): 

U)) -COs(7(t-t;)) + 

U-t r ))+ cos(7(t -*<))])}, 
where t r is the rise time defined in Section 2.2, 7 = n/t r , 

co 2 = g\k\ tanh(|/c|/i ) 

and J 7-1 is the inverse Fourier transform (see equation (3.13) below). A similar expression 
can also be derived for the velocity potential <p(x, z,t), however we do not directly need it 
in our study. 

This analytical solution will be used below in numerical simulations. It has the advan- 
tage of being simple and, thus, computationally inexpensive. However, the flat bottom 
assumption (h(x) = h = const) prevents us from using this solution beyond some small 
evolution times. The validity of this approximation has already been addressed in the 
literature [KDD07, SF09] and will be discussed at some point below. 



n =N x xN, 



77 (x, t) 



7 



i=l 



H{t-U)di{k) , 

=; [COS{u(t 

(7 2 — u 2 ) cosh(|/c|/io) 

H(t — ti — t r )[cos(oj(t 
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3.2. The weakly nonlinear (WN) model. A tsunami wave during its generation is 
usually well described by the Cauchy-Poisson problem [DD07d, KDD07, SF09]. The main 
reason for this simplification is the fact that a wave of a half meter amplitude represents only 
a tiny perturbation of over a 4000 m water column. However, the real world bathymetry 
is generally complex and may contain simultaneously various scales. For example, the 
subduction zone bathymetry represented on Figure 2 ranges from 7000 to 20 m and thus, 
nonlinear effects may be locally important. In order to take into account all realistic 
bathymetric features and study in detail the initial stages of tsunami propagation we 
describe below a new numerical model. 

We consider the physical setting and notation of Section 3.1. The governing equations 
of the classical water wave problem are the following [Lam32, Sto58, Mei94, Whi99]: 

A0 = V 2 + <9 2 > = 0, (x,z) 6 n x [-h,rj\, (3.5) 

d t r] + V0 ■ Vr] — d z <fi = 0, z — r](x, t), (3.6) 

«9 t + i|V0| 2 + i(<9 2 0) 2 + ^ = 0, z = V (x,t), (3.7) 

d t h + V0 • Vh + d z (j) = 0, z = -h(x, t), (3.8) 

with (j) the velocity potential, g the acceleration due to gravity force and V = (d x ,d y ) 
denotes the gradient operator in horizontal Cartesian coordinates. 

The assumptions of fluid incompressibility and flow irrotationality lead to the Laplace 
equation (3.5) for the velocity potential <j)(x,z,t). The main difficulty of the water wave 
problem lies on the boundary conditions. Equations (3.6) and (3.8) express the free-surface 
kinematic condition and bottom impermeability respectively, while the dynamic condition 
(3.7) expresses the free surface isobarity. 

The bathymetry h(x,t) is decomposed into the static part ho(x) (given e.g. by the 
ETOPOl database, cf. Figure 2) and the dynamic sea bed displacement ((x, t) constructed 
above in (2.1): 

h(x,t) = h (x) - ((x,t). (3.9) 

Remark 4. Recently, some weak dissipative effects have also heed included into the classical 
water wave problem (3.5) - (3.8). For more details on the vis co-potential formulation we 
refer to [DDZ08, DD07c, Dut07, Dut09b, Dut09a]. 

In the sequel we will need the unit exterior normals to the fluid domain. It is straight- 
forward to obtain the following expressions for the normals at the free surface and bottom 
respectively: 

ft ' = 7nW™' a ' = 7W [ " v ' ! '" 11 '' 

In 1968 Zakharov proposed a different formulation of the water wave problem based on 
the trace of the velocity potential at the free surface [Zak68]: 

(p(x,t) :— <p(x,T](x,t),t). (3.10) 

This variable plays the role of generalized momentum in the Hamiltonian description of 
water waves [Zak68, DB06]. The second canonical variable is the free surface elevation 77. 
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Another important ingredient is the normal velocity at the free surface v n which is 
defined as: 



V n {x,t) := y/l+WW 



{d z <j>-V<j>-Vrj)\. (3.11) 



z=r) 



The boundary conditions (3.6) and (3.7) on the free surface can be rewritten in terms of 
(p, v n and r] [CSS92, CS93, FCKG05]: 

d t r] — V v ((p) = 0, 

^ + i|V^| 2 + ^-^ T ^ T [^(^) + V^-Vr ? ] 2 =0. (3 ' 12) 

Here we introduced the Dirichlet-to- Neumann operator (D2N) V v : ip H- v n [CM85, CS93] 
which maps the velocity potential at the free surface ip to the normal velocity v n . The name 
of this operator comes from the fact that it denotes a correspondance between Dirichlet data 



ip and Neumann data \/l + | Vr/| 2 



on the free surface. We provide in Appendix B 



the complete derivation of Zakharov's formulation for the water wave problem. 

3.2.1. Numerical evaluation of the D2N operator. We saw above that the water wave prob- 
lem can be reduced to a system of two PDEs governing the evolution of the canonical 
variables rj and ip. In order to solve this system of equations we must be able to com- 
pute efficiently the quantity V v ((p). In this section we present a simple method for the 
numerical computation of the D2N operator, which is appropriate for the application of 
the linearized Euler model in the solution of problems dealing with tsunami generation. 
This approach is based on the extensive use of Fourier transforms. On the discrete level 
this transformation can be efficiently implemented with the Fast Fourier Transform (FFT) 
algorithm [CT65, FJ05]. 

The direct J 7 and inverse J 7-1 Fourier transforms in 2D are defined as follows: 

m = f0) = f f{x)e~ its dx, F-itfl = f(x) = j^y 2 J ffi)^ dk. (3.13) 

R 2 R 2 

The problem to be solved is 

+ = 0, (x,z)enx[-h,Tj\, (3.14) 

= cp, z = v , (3.15) 

y/T+ |V/i| 2 ^ = d t h, z = -h. (3.16) 

Once the function is determined, we must compute its normal derivative on the free 
surface (3.11). 

Since a tsunami wave induces a special flow regime in which the horizontal extent is 
much more important than the variations in the vertical direction, we can apply the Fourier 
transform to the Laplace equation (3.14) as if it were posed in a strip-like domain: 

d 2 d> . -*, 9 - 
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The general exact solution to this ODE can be easily computed: 

0(jfe; z) = A(k) cosh(|£|z) + B(k) swh(\k\z). (3.17) 

The two unknown functions A(k) and B(k) must be determined from the boundary con- 
ditions (3.15), (3.16). For the sake of convenience we rewrite the Neumann boundary 
condition at the bottom (3.16) in this form: 



dz 



-d t h - V0| , • Vh = f(x, t). 



(3.18) 



The right-hand side will be denoted by f(x, t), which implicitly depends on the solution <f>. 

The application of the boundary conditions leads to the following system of linear equa- 
tions: 



cosh(|/c|7/)v4(/c) + sinh(|fc|^) J B(A;) 
-\k\ smh(\k\h)A{k) + |jfe| cosh{\k\h)B{k) 
which can be easily solved: 

ip cosh(\k\h) — j 
A(k) = — 



\k\ 



(psmh{\k\h) + / 



B(k) 



cosh(| fc|?7) 



cosh(\k\H) cosh(\k\H) 

Here, H = h + 77 is the total water depth. The knowledge of these functions provides the 
velocity potential in the whole domain thanks to the general solution (3.17). 

Finally, we compute the normal velocity v n on the free surface (3.11). If we compute this 
quantity in Fourier space, the answer will be given immediately by the inverse transform 
J 7 ^ 1 . The first term of v n is readily given by the formula 



d z <j) =0\k\tanh(\k\H) + fsech(\k\H). 

Z=T) 

To compute the second term we use the following approximate expression: 

V<£\Z^Vri = fIf- 1 [ik<p] ■ J=- x [ikfj] 



(3.19) 



Remark 5. Equation (3.18) indicates that the function f(x,t n ) depends implicitly on the 
unknown solution <p(x, z,t n ) . In order to compute this apparent contradiction, we apply a 
fixed-point iteration initialized with the value of f(x, t n _i) from the previous time step: 



f 



fc+i 



•dth - 7 



v0U_ h (/*) ■ vh 



The last product is computed in the physical space: 

l (f k )-Vh}=F\F- 1 [V0 



Simple computations yield 



(f k )=ik 0sech(\k\H)-f 



/° = /(fc;t n _ 1 

(/*)] ■ V* 



k tanh(\h\H) 



\k\ 
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Our numerical experiments show that this iterative procedure is convergent and the tolerance 
e := II./ +1 ~~ ./I loo — 10~ 5 is reached after four iterations in average. 

The resulting model is only weakly nonlinear since Laplace's equation (3.14) is solved 
using the Fourier transform in a strip-like domain. Consequently, there is an implicit lin- 
earization in the solution procedure. However, the WN model contrary to the CP model 
not only takes into account some nonlinear effects but can also be efficiently applied to 
cases with realistic bathymetry. We note that this model is similar to the first order ap- 
proximation model proposed in [GN07] if in our method we further simplify all expressions 
by replacing the total water depth H by the undisturbed depth h. 



3.3. Time integration. Applying the above Fourier type spectral method to equations 
(3.12) governing the evolution of the canonical variables r] and (p leads to a system of 
ordinary differential equations, i.e. 

§ t = A(t,$), $(*„) = $o, $ = (r/,y?) T (3.20) 

In order to integrate numerically this system of ODEs we apply an integrating factor 
method analogous to the one used in [FCKG05, XG09]. This method apparently decreases 
the stiffness of the system of ODEs and therefore allows for an efficient application of 
explicit time integration schemes. We start by extracting the linear part of equations 
(3.20): 

$ t + £■ $ = (3.21) 

fo ~ — 

where C — I ^ 1 and u = y g\k\ tanh(|A;|/io) is the wave frequency corresponding 

to the wave number \k\. For a general bathymetry we choose the constant ho to be the 
mean water depth. (We note that we use the arithmetic average of values provided by the 
ETOPOl database in the region under consideration.) The term A/"($) incorporates the 
remaining nonlinear terms: 



The linear terms can be integrated exactly by the following change of variables: 
m) ._ Mt-U>)&( t ) e ^t-to) = ( cos ( w (* - *o)) -7 sin(w(t - t ))\ 

Consequently, we solve in practice the following system of ODEs: 

This simple modification allows us to take larger CFL numbers, thus improving the overall 
time stepping performance. 
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Finally, the system of ODEs is discretized by the standard fourth-order Runge-Kutta 
(RK4) scheme [HNrW09]: 

* n+ i = ^ n + \At(k l + 2k 2 + 2k 3 + k A ), 

h = B(t n ,* n ), 

k 2 = B(t n + |At,# n + ±At h), (3.22) 

k 3 = B(t n + \At,m n + \Atk 2 ), 

k A = B(t n + At,V n + Atk 3 ), 

where the subscript refers to the discrete time instance \l/ n := ^(t n ) and At is the discrete 
time step: t n+ i = t n + At. 

In the computations described below, we use a Runge-Kutta (4,5) scheme with an adap- 
tive time step control (cf. [DP80]). However it is not so fundamentally different from the 
classical RK4 scheme (3.22) described above. 

3.4. The BBM-BBM type system. When the long wave approximation is applied to 
the water wave problem (3.5) - (3.8), one obtains the well-known nonlinear shallow water 
(or Saint- Venant) equations [dSV71, Sto58, Whi99] which have been extensively used for 
tsunami simulations [Ima96, TG97, TS98, DKK08, DPD11]. If we go further in the asymp- 
totic expansions, some dispersive effects can be included and generally the resulting system 
is referred to as Boussinesq system [Bou72, BCS02, MBS03, DD07b, DMS07, DMS09]. 

In this study we use the Boussinesq system of BBM-BBM type with variable bottom 
derived in [Mit09]. See also [Per67, Cha07]. The system in dimensional variables can be 
written as: 

Vt + V • ((h + 7])u) + V • {Ah 2 [V(Vh ■ u) + Vh V ■ u\ - bh 2 Vr] t } + 

AV-(h 2 V( t ) + (t =0, (3.23) 
u t + gVri + \V\u\ 2 + Bgh [V(Vh ■ Vr?) + Vh Ar]} - dh\Au t - Bh V( tt = 0, 

where A, B, b and d are constants defined as: 

[2 2 [2 1 

A = \ , B = l-\-, b = d = -. 

V 3 3' V3' 6 

The variable u(x, t) denotes the horizontal velocity of the fluid at z = —h + a/2/3(^ + h), 
and the bathymetry variables h(x,t), h (x), C(^;^) are defined in Section 2. 

We integrate numerically the system (3.23) by using the standard Galerkin/finite el- 
ement method with PI elements for the spatial discretization coupled with an explicit, 
second-order Runge-Kutta method for the temporal discretization (so-called improved Eu- 
ler scheme) [HNrW09]. A proof that the semidiscrete system is not stiff and thus that the 
specific RK method is sufficient can be found in [DMS10]. 

In order to obtain a well-posed problem, we impose homogeneous Dirichlet boundary 
conditions which absorb partially the wave while reflecting only small amplitude oscilla- 
tory waves. Moreover, the specific numerical method appears to converge with optimal 
rate in the L 2 and L°° norms whether we consider structured or unstructured grids. This 
is contrary to the analogous initial boundary value problems with zero Dirichlet boundary 
conditions on u for the Peregrine system [Per67] where the analogous numerical method 
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converges with suboptimal orders on structured and unstructured grids. For more informa- 
tion on the properties and the implementation of the numerical method for a BBM-BBM 
type system we refer to [DMS07, Mit09]. 



4. Numerical results 

In this section we compare the propagation of a solitary wave when it is used as an 
initial condition in both the CP and WN models. Moreover, we study the generation and 
the initial stages of the propagation of the tsunami wave of the July 17, 2006 event. We 
also present a comparison between the WN, CP and Boussinesq models. 



4.1. Solitary wave propagation. Before performing the Java 2006 tsunami generation 
simulations, we study the propagation of a solitary wave solution to the full water wave 
problem using the WN and CP models. The initial condition is a solitary wave, computed 
by using the method presented by Tanaka [Tan86]. 

Consider the two-dimensional water wave problem in a channel of uniform depth ho = 
const. Since we look for travelling wave solutions, the flow field can be reduced to the steady 
state by choosing a frame of reference moving with the wave speed c. The introduction of 
dimensionless variables leads to a single scaling parameter, the Froude number Fr defined 

as Fr := — =. Hereafter, the governing equations are considered in dimensionless form. 

Vgh 

The complex velocity potential is classically introduced as w = <fi + iip, where ip is the 
stream function. We choose = at the crest and if = at the bottom. The fluid 
region is then mapped onto the strip < if < 1, — oo < < oo on the plane w with 

ip — 1 corresponding to the free surface. We introduce the quantity Q = log — — = r — i9, 

dz 

where 9 is the angle between the velocity vector and horizontal axis Ox. The real part r is 
expressed in terms of the velocity magnitude q as r = log q. The boundary conditions to 
be satisfied are the dynamic condition on the free surface and the bottom impermeability 
which are expressed as 

dq 3 3 

— — = ~ sin 9, on ip — 1 and = 0, on ip = 0. (4.1) 

dcp Fr 

Consequently, the problem is now transformed into the determination of the complex 
function Q, analytic with respect to w within the region of the unit strip < if < 1, 
decaying at infinity and satisfying the boundary conditions (4.1). By applying Cauchy's 
integral theorem, one can find the following integral equation on the free surface ip = 1: 

00 00 00 

2 f m jfi^ ; t [tWL 



7T J (f — 0) 2 +4 71 J (if — 0) 2 +4 7T J If 

—00 —00 —00 



where r(0) and 8((f>) denote the traces of the corresponding functions on the free surface 
if = 1. 
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Figure 6. Solitary wave solutions of various amplitudes for the full water 
wave problem. Both x and rj have been non-dimensionalized by the depth 
ho. 



The integral equation is solved iteratively. The convergence is tested with respect to 
the Froude number. Several solitary wave solutions computed in this way are plotted on 
Figure 6 for illustrative purposes. 

In order to illustrate the advantages of the proposed WN model over the classical CP 
solution, we let a solitary wave with amplitude A/h = 0.1 propagate up to T = 80 (in this 
section we use dimensionless quantities and time T is non-dimensionalized by ^/g/h ). 

We recall that the classical CP solution of (3.1)-(3.4) corresponding to the initial free 
surface height r]\ t=0 = rj (x) and the velocity potential distribution at the free surface 
f\t=o = fo{ x )i takes the following form: 

i](x, t) = J 7 " 1 |^o(^) cos(a;t) H — <fio(k) sin(ut) j , 

<f)(x, z,t) = J 7-1 1 {(po{k) cos(ut) — —f)o(k) sin(u;t)) (cosh(|£|2:) + tanh(|fc|/i) sinh(|fc|z)) j, 

where f)o(k) = J-"{?7o(x)} and (po(k) = ^{^^{x)} are the Fourier transforms of the initial 
conditions. 

The solution profiles of both models are presented in Figures 7 (a)-(e). One observes 
that the WN model preserves quite well the shape of the solitary wave while shedding a 
small dispersive tail behind. The CP solution gradually transforms the initial wave into a 
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dispersive tail according to the linear nature of equations (3.1)-(3.4). In Figure 7 (f) we 
present the normalized amplitude error defined as: 

| max{r)(x, £)} — A/h \ 



where ma,x{r](x,t)} denotes the discrete maximum of the numerical solution and Aj ho = 

0.1 is the exact solitary wave amplitude. In both computations a uniform grid of 512 nodes 
is used. Here, again, we notice a better performance of the WN solver compared to that 
of the CP solution. This specific experiment shows that the WN model is a better model 
compared to the CP solution when nonlinear effects must be included for the study of 
tsunami generation and propagation. 



4.2. The July 17, 2006 tsunami generation simulation. The main purpose of this 
study is to present a novel methodology for tsunami generation problems. This approach 
is illustrated on the example of the July 17, 2006 Java tsunami since this event is not 
completely understood yet and there is an available finite fault solution for the presumed 
generating underwater earthquake. 

In this section we show a practical application of the WN method for water waves 
generated by a moving bottom. Namely, we exploit the bottom motion (2.1) constructed 
in Section 2.2. The corresponding hydrodynamic problem is solved by the three methods 
discussed above: the linearized water wave problem (CP), BBM-BBM system and the 
novel WN model. 

The solution given by the WN model and the exact solutions to the linearized Euler 
equations (3.1) - (3.4) are computed on a uniform grid of 512 x 512 points. The time 
step At is chosen adaptively according to the RK(4,5) method proposed in [DP80]. The 
BBM-BBM system is solved on a triangular unstructured grid of 86276 elements. The time 
integration is performed with the classical RK2 scheme [HNrW09] with time step At = 0.5 
s. 

Several snapshots of the free surface elevation computed with the WN model are shown 
in Figures 8 (a) - (f). Analogous contour plots of the solutions of the CP and BBM-BBM 
models are almost identical and differences cannot be observed within graphical accuracy. 
Therefore, they are not presented here. The parameters of the bottom motion, bathymetry 
and computational domain geometry were explained in Section 2. 

In this computation, we see a complex process of simultaneous wave evolution together 
with rupture propagation during approximately 210 s. Namely, the free surface deformed 
by the rupture of the first subfaults evolves while the rupture continues to propagate along 
the fault. This kind of fluid/moving bottom interaction cannot be described in the static 
generation framework, cf. Figure 10. 

In order to compare the three models described above we put eight numerical wave 
gauges at the following locations: six close to the source ((a) (107.2°, —9.388°), (b) (107.4°, 
-9.205°), (c) (107.6°, -9.648°), (d) (107.7°, -9.411°), (e) (108.3°, -10.02°), (f) (108.2°, 
-9.75°)) and two further away from the source area ((g) (108°, -10.5°), (h) (108°, -9°)). 
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Figure 7. Propagation of a solitary wave with the weakly nonlinear method 
(solid line) and Cauchy-Poisson solution (dashed line). The solitary wave 
amplitude is A/ho = 0.1. Space has been scaled by ho and time by ygjho- 



The locations of the wave gauges are represented by the symbol o on Figure 9 along with 
the static sea bed displacement. 

The eight wave gauge records are presented in Figures 10 (a)-(h). In order to show 
the importance of the dynamics of the rupture process, the records obtained from the 
static approach are also included. The overall agreement among the three dynamic models 
appears to be satisfactory. We underline that the CP solution is very close to the other 
solutions despite the fact that the bathymetric features are neglected. We also note that 
the specific BBM-BBM type system underestimates by a small amount the maximum wave 
amplitude compared to the WN model. Further numerical tests showed some sensitivity of 
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Figure 8. Snapshots of the free surface elevation computed with the weakly 
nonlinear (WN) model. Water waves are generated by dynamic co-seismic 
bottom displacements (2.1) reconstructed using the corresponding finite fault 
solution [Ji06]. 
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Figure 9. Location of the eight numerical wave gauges (indicated by the 
symbol o) superposed with the static co-seismic bottom displacement. 

the BBM-BBM solution to the bottom motion scenario [DD07d]. Namely, we can report, 
for example, that the exponential scenario led to a slighty larger wave amplitude compared 
to the other models. As expected, the static approach exhibits differences both in the 
shape and in the arrival time of the waves. Further away from the source area, the CP 
solution continues to be accurate. This is due to the fact that nonlinearity is not important 
during the propagation stage of such small amplitude waves. 

5. Conclusions and perspectives 

In the present work we considered an important issue in the modeling of tsunami gener- 
ation. Namely, a new method for the construction of dynamic co-seismic sea bed displace- 
ments was proposed. This method basically relies on two main ingredients: 

• the finite fault solution [BLM00, JWH02] gives the slip distribution along the fault 

• dynamic sea bed deformation scenarios [Ham73, DDK06, DD07d] allow us to take 
into account available information of the rupture dynamics 

To our knowledge, this reconstruction of the bottom motion is new. All developments 
presented in this paper are illustrated on the example of the July, 17 2006 Java event. 

Along with the bottom motion construction, we discussed three models to solve approx- 
imately the corresponding hydrodynamic problem and compute the induced free surface 
motions. The July 17, 2006 tsunami generation case was computed with three different 
models and a comparison was performed. We obtained a surprisingly good agreement 
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Figure 10. Free surface elevation computed numerically with four models 
at eight wave gauges located approximately at the local extrema of the static 
bottom displacement. The elevation (vertical axis) is expressed in meters, 
while the time (horizontal axis) is in seconds. The models 1, 2 and 4 use 
the dynamic finite-fault rupture (weakly nonlinear model, linearized Euler 
equations, BBM-BBM model). The third model uses the static approach 
(weakly nonlinear model). 



22 



D. DUTYKH, D. MITSOTAKIS, X. GARDEIL, AND F. DIAS 



between the CP solution and the solutions of the other two models. Recall that in the 
latter the bottom is assumed to be flat. Discrepancies will appear later in time since the 
bathymetry plays a crucial role in the tsunami propagation. 

Taking into account the simplicity and the relatively good accuracy of the new WN 
approximation to the full water wave problem with time dependent variable bottom, we 
suggest its use for the computation of the initial stages (~ 300 s) of the life of a tsunami. 
The propagation and runup can be computed afterwards by other sophisticated tools 
[TG97, IYO06, SBT + 07, DPD11], some of them being already integrated into tsunami 
warning systems [TGB + 05, WL08]. 

However we point out that extreme runup values measured after the July, 17 Java 2006 
event [FKM + 07] deserve additional numerical studies. 

Appendix A. Finite fault parameters 



Table 3: Subfault parameters given by the finite fault inver- 
sion [Ji06] . 
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107 


58636 


12 


74141 


112.11020 


99.51801 


-9 


44385 


107 


45708 


12 


74141 


60.23628 


97.77266 


-9 


40006 


107 


32778 


12 


74141 


126.96870 


80.27277 


-9 


35628 


107 


19849 


12 


74141 


63.39000 


65.00801 


-9 


31249 


107 


06921 


12 


74141 


0.52621 


94.79313 


-9 


26871 


106 


93991 


12 


74141 


1.52171 


66.78681 


-9 


22492 


106 


81062 


12 


74141 


10.96743 


81.94861 


-9 


18114 


106 


68134 


12 


74141 


2.38062 


123.04830 


-9 


96479 


109 


29915 


14 


71768 


22.40949 


123.90350 


-9 


92100 


109 


16986 


14 


71768 


48.62879 


115.45630 


-9 


87722 


109 


04057 


14 


71768 


5.99559 


83.81007 


-9 


83343 


108 


91128 


14 


71768 


7.22945 


123.80940 


-9 


78965 


108 


78199 


14 


71768 


0.10031 


93.40998 


-9 


74586 


108 


65269 


14 


71768 


0.36991 


69.37087 


-9 


70208 


108 


52341 


14 


71768 


104.18760 


123.83230 


-9 


65829 


108 


39411 


14 


71768 


46.12533 


95.97049 


-9 


61451 


108 


26482 


14 


71768 


0.28679 


89.56866 


-9 


57072 


108 


13554 


14 


71768 


2.06597 


80.14312 


-9 


52694 


108 


00624 


14 


71768 


30.55070 


66.23147 


-9 


48315 


107 


87695 


14 


71768 


73.72994 


87.91253 


-9 


43937 


107 


74767 


14 


71768 


112.90700 


92.28181 


-9 


39558 


107 


61837 


14 


71768 


74.73608 


86.51558 


-9 


35180 


107 


48908 


14 


71768 


121.73820 


64.68654 


-9 


30801 


107 


35979 


14 


71768 


231.20940 


65.50779 


-9 


26423 


107 


23050 


14 


71768 


96.55727 


87.01543 


-9 


22044 


107 


10121 


14 


71768 


28.29534 


122.55670 


-9 


17666 


106 


97192 


14 


71768 


0.84110 


70.21989 


-9 


13287 


106 


84263 


14 


71768 


7.99213 


87.51706 


-9 


08909 


106 


71334 


14 


71768 


1.33281 


96.33266 


-9 


87274 


109 


33115 


16 


69394 


43.31154 


121.79150 


-9 


82896 


109 


20187 


16 


69394 


87.17052 


124.49750 


-9 


78517 


109 


07257 


16 


69394 


61.47630 


87.10537 


-9 


74139 


108 


94328 


16 


69394 


31.53286 


70.58137 


-9 


69760 


108 


81400 


16 


69394 


0.70628 


65.17896 


-9 


65382 


108 


68470 


16 


69394 


5.74160 


87.70702 


-9 


61003 


108 


55541 


16 


69394 


93.47714 


107.32000 


-9 


56625 


108 


42612 


16 


69394 


93.55753 


85.39201 


-9 


52246 


108 


29683 


16 


69394 


47.25525 


74.24297 



Continued on the next page 
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Latitude, ° 


Longitude, ° 


Depth, km 


Slip, cm 


Rake, ° 


-9.47868 


108.16754 


16 


69394 


24.65230 


124.20110 


-9.43489 


108.03825 


16 


69394 


35.63115 


71.78733 


-9.39111 


107.90896 


16 


69394 


25.11757 


75.27779 


-9.34732 


107.77967 


16 


69394 


68.15302 


107.42980 


-9.30354 


107.65038 


16 


69394 


24.66007 


112.77880 


-9.25975 


107.52109 


16 


69394 


0.50688 


79.86887 


-9.21597 


107.39180 


16 


69394 


119.92850 


75.03103 


-9.17218 


107.26250 


16 


69394 


77.08335 


110.83160 


-9.12840 


107.13322 


16 


69394 


31.65430 


123.83060 


-9.08461 


107.00393 


16 


69394 


11.42768 


66.47282 


-9.04083 


106.87463 


16 


69394 


33.80650 


115.65650 


-8.99704 


106.74535 


16 


69394 


39.47481 


65.15574 


-9.78069 


109.36316 


18 


67021 


35.42621 


111.95830 


-9.73691 


109.23387 


18 


67021 


103.05030 


124.62650 


-9.69312 


109.10458 


18 


67021 


101.38220 


122.70620 


-9.64934 


108.97529 


18 


67021 


76.76701 


68.20042 


-9.60556 


108.84600 


18 


67021 


10.71945 


77.79713 


-9.56177 


108.71671 


18 


67021 


1.32449 


100.72950 


-9.51799 


108.58742 


18 


67021 


37.46857 


124.59330 


-9.47420 


108.45813 


18 


67021 


118.99580 


100.38000 


-9.43042 


108.32883 


18 


67021 


79.62616 


91.56905 


-9.38663 


108.19955 


18 


67021 


97.61735 


109.86430 


-9.34285 


108.07026 


18 


67021 


87.67753 


87.57239 


-9.29906 


107.94096 


18 


67021 


15.14859 


64.75201 


-9.25528 


107.81168 


18 


67021 


82.60960 


71.66805 


-9.21149 


107.68239 


18 


67021 


66.06397 


98.55843 


-9.16771 


107.55309 


18 


67021 


0.43085 


67.81042 


-9.12392 


107.42381 


18 


67021 


35.30429 


124.04570 


-9.08014 


107.29452 


18 


67021 


59.17323 


124.55130 


-9.03635 


107.16522 


18 


67021 


15.23214 


66.82615 


-8.99257 


107.03593 


18 


67021 


28.10358 


76.08198 


-8.94878 


106.90664 


18 


67021 


48.09923 


124.24450 


-8.90500 


106.77735 


18 


67021 


42.38682 


124.42850 



Appendix B. Zakharov's formulation of the water wave problem 

In this appendix we recast the governing equations (3.5) - (3.8) of the water wave 
problem in a more compact and mathematically more convenient form [Zak68, CS93]. 

Using the definition of the normal velocity (3.11), it is straightforward to rewrite the 
kinematic free surface condition (3.6): 

d t rj - V v ((p) = 0, 

where (p is the trace of the velocity potential at the free surface (3.10). 



V#U = V*. - V, «U = "»"' , (B.5) 
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The time derivative and the horizontal gradient of the velocity potential trace on the 
free surface can be computed: 

d t ip = d t <\> + d tV d z <f)\ z=ri = dt<f> + V v (ip) d z <j)\ z=v , (B.l) 

and similarly one can compute the horizontal gradient: 

= V0| z=r? + Vr] d z (f)\ z=fl . (B.2) 

In order to close the system, we have to express all derivatives of the potential computed 
at the free surface, in terms of ip, r\ and T> v (tp). 

From the definition of the normal velocity (3.11) and the D2N operator one readily 
obtains: 

S/4>\ z=ri -V V =d z 4>\ z ^-V T ,(cp). (B.3) 
Substituting the last identity into (B.2) multiplied by V77, leads to the following expression: 

1 + |Vr/| 2 

Now we have all elements to find the horizontal derivatives of the velocity potential: 

;i + |Vr/| 2 )Vy? - Vr, ((p)Vy - (V<p ■ Vr))Vr) 

In order to rewrite Bernoulli condition (3.7) in new variables, we make the following 
observation (using (B.2) and (B.3)): 

||V0| 2 + \{d z ct>) 2 = |V0 • + \d z ct> d z <p = 

= ±V0 • (Vip - d^Vri) + id t (f>(V v (ip) + V0 • V77) = ±V0 • + lV v (ip)d z( j), z = V . 

Taking into account this observation and expression (B.l) for the time derivative of (p, the 
dynamic condition takes this equivalent form: 

d t if + grj + |V0 • S7ip - \U v {ip)d z (j) = 0, z = rj. 

After substituting expressions (B.4), (B.5) into the last equation and summarizing all the 
developments made above, we get the following set of dynamic equations equivalent to the 
complete water wave problem (3.5) - (3.8): 

d t 7] - V v (ip) = 0, 

d t ip + || Vyf +g V - ^TTWW) l V M + ^ - Vr/] 2 = 0. 

Appendix C. Relations between elastic constants 

In the classical elasticity theory, coefficients in Lame equations (governing the displace- 
ments field in an elastic solid), can be expressed in terms of various sets of physical pa- 
rameters [Lov44, SS46]. The purpose of this Appendix is to recall some relations between 
them. 
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Lame coefficients A and /i can be denned in terms of the Young's modulus E (having the 
dimension of the pressure [Pa]) and Poisson's ratio v (dimensionless coefficient < v < 
1/2): 

Ev E 



A 



;i + 2i/) 



and inversely: 



(3A + 2/i)/i 
hi = — , v 



2(l + i/) 
A 



A + /z ' 2(A + /i)' 

The celerities of P and 5 waves have the following expressions in terms of Lame coefficients: 



f A + 2/i 

j c s — 

P 

where p is the density of elastic medium. These relations yield 

V- = P c li * = P c l ~ 2 P- 
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